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ABSTRACT 

The  use  of  advanced  techniques  can  greatly  improve  the 
effectiveness  of  Monte  Carlo  simulation  calculations.  As  a 
demonstration  model,  the  Navy's  Antisubmarine  Warfare  Air 
Engagement  Model,  APAIR,  which  simulates  a  single  aircraft 
hunting  and  destroying  a  submarine,  was  selected.  Possible  improve¬ 
ments  in  random  number  generation  are  presented;  however,  the 
study  centers  on  implementation  of  variance  reduction  techniques. 

Two  test  cases,  typical  of  APAIR  implications,  were  chosen. 

Examples  illustrating  the  use  of  the  statistical  estimation,  expected 
value,  systematic  sampling,  antithetic  sampling,  correlated  sampling, 
history  reanalysis,  and  importance  sampling  were  run.  A  comparison 
of  variances  with  the  unmodified  APAIR  showed  the  effectiveness  of 
variance  reduction  techniques.  Efficiencies,  equivalent  to  reduced 
running  time  to  obtain  the  same  variance,  of  nearly  a  factor  of  20 
were  obtained  in  specific  cases.  It  was  felt  that  throung  experience 
and  careful  effort  overall  improvements  of  a  factor  of  10  could  be 
expected. 
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EXECUTIVE  SUMMARY 


Contemporary  design  and  construction  of  large  scale  Monte  Carlo 
systems  analysis  programs  seldom  consider  incorporation  of  the  efficient 
simulation  techniques  that  have  been  developed,  tested  and  proven  successful 
in  various  technical  disciplines.  Most  notable  among  these  are  random 
number  generation  and  variance  reduction  schemes  that  have  been  routinely 
used  in  radiation  transport  to  provide  vast  improvements  in  Monte  Carlo 
simulations. 


One  objective  of  a  project  sponsored  at  Science  Applications,  Inc. 
(SAI)  by  Code  462  of  the  Office  of  Naval  Research  was  to  develop  these  tech¬ 
niques  to  the  point  where  they  would  be  generally  applicable.  Basically,  this 
involved  development  of  improved  techniques  for  selecting  probability  distri¬ 
butions,  schemes  for  generation  of  random  numbers  and  procedures  for 
variance  reduction  for  general  application  in  Monte  Carlo  simulation.  The 
results  of  these  developments  are  presented  in  Refs.  1,  2,  and  3. 


Another  objective  of  the  project  was  to  demonstrate  their  applica¬ 
bility  to  a  large  scale  Navy  simulation  program. 

It  is  the  purpose  of  this  document  to  summarize  the  results  of  the 

demonstration  effort  for  the  Antisubmarine  Warfare  Air  Engagement  Model 
(4  5  6) 

A  PAIR.  ’  ’  The  efficiency  estimates  showed  that  significant  improvements 
in  APAIR  running  times  (about  a  factor  of  20  in  some  cases)  can  be  achieved 
using  the  improved  simulation  techniques  investigated  during  this  project. 

It  was  generally  felt  that  by  judicial  selection  of  variance  reduction  schemes, 
that  running  times  for  the  same  accuracy  can  be  reduced  by  a  factor  of  10 
for  many  APAIR  problems. 
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1.  INTRODUCTION 


The  research  performed  in  a  project  sponsored  by  the  Office  of 
Naval  Research  (Code  462)  at  Science  Application,  Inc.  (SAI)  over  die  past 
year  was  directed  toward  achieving  the  following  objectives: 


•  Develop  improved  probability  selection  techniques. 

•  Develop  improved  random  number  generation  procedures  for 
selected  probability  distributions . 

•  Improve  variance  reduction  technology. 

•  Demonstrate  the  application  of  improved  simulation  techniques 
to  the  ASW  air  engagement  simulation  model  APAER.  °) 

The  results  of  the  effort  performed  in  the  first  three  areas  above  are  docu¬ 
mented  in  Refs.  1,  2  and  3  respectively.  It  is  the  purpose  of  this  document 
to  summarize  the  results  of  the  effort  that  involved  demonstration  of  improved 
Monte  Carlo  Simulation  techniques  for  the  APAIR  ASW  engagement  model. 

The  APAIR  model,  selected  for  the  study  here,  is  a  Monte  Carlo 
simulation  of  the  full  engagement  between  one  aircraft  and  one  submarine. 

The  general  version  (APAIR  2. 6)  simulates  the  detection,  localization  and 
attack  phases  of  airborne  ASW  missions.  The  model  is  particularly  suited 
for  this  studv  since  it  is  a  relatively  long  running  program,  and  because  it 
is  in  rather  wide  use,  improvements  in  the  running  dme  would  be  most 
welcome.  It  is  emphasized,  however,  that  the  objective  was  not  to  critique 
the  APAIR  model  but  rather  to  evaluate  the  effectiveness  of  efficient  tech¬ 
niques  which  are  not  commonly  used  in  Monte  Carlo  simulations.  It  should 
also  be  mentioned  that  not  only  are  the  techniques  studied  here  applicable  to 
the  APAIR  program  but  also  to  other  ASW  simulation  models  such  as 
APSURV,  APSUB  and  APSURF. 
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The  version  of  APAIR  used  in  the  study  was  provided  to  SAI  by  Code 
141,  Naval  Undersea  R/D  Center  (NUC)  in  San  Diego.  In  addition  to  the 
program,  Code  141  also  provided  SAI  with  a  ’'representative"  set  of  data 
for  demonstration  purposes  only.  No  attempt  was  made  to  interpret  any 
of  the  results  beyond  that  necessary  to  evaluate  the  improved  simulation 
efficiency.  However,  the  results  are  presented  in  sufficient  detail  that 
this  document  should  be  generally  useful  to  the  analyst  interested  in  impro\  - 
ing  the  random  number  generation  schemes  and  incorporation  of  variance 
reduction  in  APAIR. 

Of  the  three  areas  considered  for  improvement  in  APAIR  (probab¬ 
ility  distribution  selection,  random  number  generation  and  variance 
reduction),  variance  reduction  proved  to  be  the  most  successful.  Improve¬ 
ment  in  the  probability  distribution  modeling  was  restricted  due  to  the  lack 
of  sufficient  data  on  input  parameters  (such  as  aircraft  navigation  errors, 
sonobuoy  fixes  and  drop  errors,  aircraft  weapons  effects,  tactics,  etc.). 

In  the  area  of  random  number  generation,  several  modifications 
were  identified  which  would  not  only  improve  the  generation  of  random  num¬ 
bers  but  also  increase  the  ease  with  which  random  numbers  could  be  generated 
on  various  computers.  The  latter  was  accomplished  by  developing  a  random 
number  generator  that  was  machine  independent.  This  is  particularly  impor¬ 
tant  when  comparison  between  results  from  different  machines  or  facilities 
are  desired  or  when  the  program  is  to  be  made  operational  at  new  facilities. 

The  third  area  involved  the  application  of  variance  reduction  techniques. 
This  proved  to  be  extremely  fruitful.  The  variance  reduction  techniques  used 
here  were  developed  primarily  for  application  to  radiation  transport  problems 
and  have  not  been  widely  used.  However,  it  was  shown  conclusively  here 
that  a  much  broader  application  is  warranted.  For  example,  in  some  cases 
improvements  in  efficiency  (i.  e. ,  the  running  time  required  to  achieve  the 
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same  variance)  was  a  factor  of  almost  20.  An  overall  improvement  of  a 
factor  of  10  could  generally  be  expected  if  the  effort  to  understand  and  care¬ 
fully  apply  the  various  techniques  is  expended.  A  summary  of  the  improve¬ 
ment  in  efficiency  in  terms  of  running  time  reduction  is  shown  in  Table  X.  1 
for  the  various  variance  reduction  techniques  applied  to  APAIR. 

In  the  following  section  of  the  report  the  improvements  in  random 
number  generation  are  described.  The  use  of  variance  reduction  in  APAIR 
will  be  described  in  detail  in  Section  3. 

Since  the  APAIR  study  was  for  demonstration  purposes  only,  it  should 
be  mentioned  that  only  expedient  modifications  were  made  to  the  program. 
Therefore,  the  APAIR  program  at  SAI  containing  the  improvements  discussed 
:.s  not  considered  to  be  a  version  that  would  be  generally  useful.  However, 
it  has  been  proposed  that  these  modifications  be  made  to  APAIR  to  provide 
for  user  flexibility  and  improved  efficiency  in  the  future. 


A  Summary  of  Efficiency  Improvements  in  APAIR  Application  Using  Variance 

Reduction  Techniques 


2.  RANDOM  NUMBER  GENERATION  IMPROVEMENTS  IN  APAIR 


In.  this  section,  the  potential  improvements  in  the  APAIR  random  num¬ 
ber  generators  are  discussed.  It  is  not  anticipated  that  replacement  of  the 
random  number  gem  rators  currently  u;ed  in  APAIR  would  significantly  re¬ 
duce  the  total  program  running  time  in  most  cases,  since  the  fraction  of 
time  spent  in  generating  random  numbers  in  a  typical  application  is  not 
large.  However,  in  some  cases,  the  approaches  recommended  here  will 
be  faster,  more  accurate  and  more  convenient. 

Possibly  the  most  significant  improvement  for  APAIR  random  number 
generation  would  be  in  the  us*  of  a  machine  independent  uniform  random  num¬ 
ber  generator.  This  will  automatically  permit  identical  sequences  of  random 
numbers  to  be  generated  on  different  computers.  The  main  advantage  in  this, 
of  course,  is  that  identical  results  can  be  obtained  on  different  machines  and 
at  different  installations  for  comparison  purposes.  Furthermore,  since  the 
existing  random  number  generator  now  used  in  APAIR  is  machine  dependent, 
transferring  the  code  to  other  machines  can  cause  considerable  difficulty. 

Use  of  the  machine  independent  version  described  here  would  eliminate  this 
problem.  The  use  of  the  machine  independent  uniform  random  number  genera¬ 
tor  and  improvements  in  the  normal  distribution  will  be  discussed  below. 

2. 1  UNIFORM  RANDOM  NUMBER  GENERATOR 

The  uniform  random  number  generator  in  the  current  version  of  APAIR 
obtained  by  SAI  from  NUC  and  designed  for  use  on  the  Univac  1108  is  a  mixed 
congruential  generator  employing  the  algorithm 

X  j  =  27095269935- Xn-f  2049  (mod  235) 

where  and  are  successive  random  numbers.  This  generator  has 

not  been  subjected  to  the  Coveyou-MacPherson  analysis  (Ref.  9),  the  most 
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exacting  test  for  random  number  generators  currently  known .  Therefore,  its 
validity  as  a  random  number  generator  is  not  proven.  There  are,  however, 
no  a  priori  reasons  to  suspect  that  this  generator  is  faulty. 

This  generator  is,  however,  not  machine  independent.  It  will  work 
only  on  the  Univac  1108  and  other  machines  with  a  36-bit  word  length  and  a 
similar  negative  integer  representation.  It  will  not  work,  for  example,  on 
an  IBM  360  computer.  A  better  choice  of  a  basic  random  number  generator 
would  be  the  machine  independent  generator,  MIRAN,  described  in  Appendix 
A*  This  would  allow  greater  flexibility  in  exchange  of  program  and  problems 
between  different  computer  facilities  and  in  addition  uses  an  algorithm  of 
proven  validity.  Random  number  generation  occupies  such  a  small  portion 
of  the  overall  APAIR  computing  time  that  the  loss  of  efficiency  entailed  in 
using  MIRAN  would  not  cause  an  increase  in  the  total  running  time  of  APAIR 
problems. 

2. 2  NORMAL  RANDOM  NUMBER  GENERATOR 

The  normal  distribution  is  used  in  APAIR  to  generate  random  varia¬ 
bles  such  as  aircraft  navigation  errors  (inertial,  tactical,  navigation,  reset, 
areas  of  uncertainty  for  sonobuoy  fixes,  etc. ).  The  APAIR  procedure  currently 
used  to  obtain  a  normal  random  variable  from  the  distribution: 


f(x) 


(2.1) 


is  to  generate  an  approximate  normal  random  deviate  using 


' 1  r-  - 


Where  Ru^, . . .  >*^2  is  a  set  of  12  independent  random  values  from  the 
uniform  distribution  U(0, 1).  This  procedure  takes  about  105  microseconds 
on  a  UNIVAC  1108  using  assembly  language. 

The  result  used  in  (2. 2)  is  based  on  the  central  limit  theorem^^ 
and  is,  therefore,  approximate.  For  example,  the  range  of  X  using  (2. 2) 
is  limited  to 

u-6a<X<fi  +  6o-  (2. 3) 

which  is  probably  adequate  in  most  situations.  Difficulty  could  arise  when 
very  small  probability  events  (i.  e. ,  outside  the  range  of  X  as  given  by 
(2. 3))  are  important. 

A  better  method  for  generating  random  numbers  from  the  normal 
distribution  which  is  both  exact  (within  machine  accuracy)  and  requires  only 
30  microseconds  is  a  technique  developed  by  Marsaglia^  and  documented 
in  Ref.  2.  A  routine  using  this  approach,  is  shown  below  along  with  the 
corresponding  flow  diagram  in  Fig.  2.1. 
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Normal  Distribut: 


Sample  Routine 


FUNCTION  RANORM  (DUMMY) 

R  -  RANUMB(R) 

IF  (R  GT.O.  8638)  GO  TO  10 

RANORM  =  2.  *(RANUMB(X)  +  RANUMB(Y)  +  RANUMB(Z)  -  1.  5) 
RETURN 

10  IF  (R  GT.O.  9745)  GO  TO  20 

RANORM  =  1.  5*(RANUMB(X)  +  RANUMB(Y)  -  1, 0) 

RETURN 

20  IF  (R  GT.O.  997302039)  GO  TO  100 
25  X  =  6.  *RANUMB(X)  -  3.  0 
Y  =  0.  358*RANUMB(X1 
XSQ  =  X*X 

GX  =  17. 49731196*EXP(-XSQ*.  5) 

AX  -  ABS(X) 

IF  (AX.  GT.  1. 1)  GO  TO  30 

IF  (Y.GT.(GX- 17. 44392294  +  4. 73570328+XSQ  +  2. 157987544*AX)) 
GO  TO  25 
RANORM  =  X 
RETURN 

30  AX3  =  2. 367985163*(3~AX)**2 

IF  (AX.  GT.  1.  5)  GO  TO  40 

IF  (Y.  GT.  (GX-AX3-2. 157987544*(1.  5-AX)))  GO  TO  25 

RANORM  =  X 

RETURN 

40  IF  (Y.  GT.  (GX-AX3))  GO  TO  25 
RANORM  =  X 
RETURN 

100  X  =  SQRT  (9-2*(ALOG(RANUMB(X))) 

IF  (RANUMB(X).  GT.  3/X)  GO  TO  100 

IF  (RANUMB(X).  GT.  0.  5)  X  =  -X 

RANORM  =  X 

RETURN 

END 

As  written  above,  the  routine  generates  a  normal  standard  deviate 
(i.  e. ,  a  random  number  from  (2. 1)  with  cr  =  1  and  /x  =  0)  defined  as  N(0, 1). 

It  is  then  left  up  to  the  calling  program  to  multiply  by  the  standard  deviation 
and  add  the  mean  if  a  generalized  normal  deviate  is  required.  That  is,  for 
a  distribution  with  mean  n  and  variance  cr%}  (i.  e. ,  Eq.  2. 1)  the  correct 
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random  number  would  be  crN(0, 1)  +  p  ,  where  N(0, 1)  is  a  random  number 

2 

from  a  distribution  with  p  =  0  and  a  =  1. 

Therefore,  if  the  algorithm  prescribed  above  is  used  in  APAIR ,  the 
running  time  for  generation  of  random  numbers  from  a  normal  distribution 
would  improve  almost  a  factor  of  4  and  it  also  will  provide  exact  answers. 


■jm*  K  ■  i- Hi  i  ^mwjrw,  ujf^j^m  i  if  Pll  i  *  **^  *  r* 


3.  VARIANCE  REDUCTION  IMPROVEMENTS  IN  APAIR 


The  type  of  problems  solved  by  APAIR  provide  an  excellent  opportunity 
to  apply  a  wide  variety  of  variance  reduction  techniques  which  can  substan¬ 
tially  improve  APAIR  efficiency.  Several  of  these  were  tried  on  APAIR  with 
generally  excellent  results.  In  some  applications  the  gain  in  efficiency  was 
estimated  to  be  almost  a  factor  of  20.  It  appears  that  a  factor  of  10  improve¬ 
ment  can  readily  be  accomplished  in  many  APAIR  applications  if  the  time 
and  effort  is  expended  to  understand  and  implement  the  appropriate  techniques. 

It  is  unquestionably  to  the  benefit  of  the  analyst  to  follow  such  an  approach 
when  long  running  times  are  of  concern. 

Resources  available  to  perform  the  study  did  not  permit  application  of 
all  the  possible  variance  reduction  schemes  available.  Therefore,  a  selected 
number  of  techniques  representing  a  broad  range  of  characteristics  were  applied 
to  various  types  of  ASW  problems.  These  included: 

•  Statistical  Estimation  using  an  expected  value  of  kill  at  each 
torpedo  drop  rather  than  scoring  actual  kills  based  on  Monte 
Carlo  simulation.  Two  cases  were  studied  corresponding  to 
different  types  of  Magnetic  Anomaly  Detectors  (MAD). 

•  Expected  Value  replacing  actual  submarine  kill  by  a  reduction 
in  the  "survival  value"  of  the  submarine  for  computing  percent 
kill  as  a  function  of  which  torpedo  caused  the  kill  (first,  second, 
etc. ).  This  was  also  performed  for  two  types  of  MAD  gear. 

•  Systematic  Sampling  selecting  initial  starting  coordinates  for 
the  submarine  location  using  two  types  of  MAD  gear.  Two  types 
of  systematic  sampling  were  considered. 

e  Antithetic  Variates  used  to  select  initial  starting  locations  of  the 
submarine. 

•  Correlated  Sampling  using  random  number  control  to  generate 
identical  histories  to  a  point  where  differences  in  two  types  of 
MAD  gear  lead  to  a  divergence  of  histories.  This  was  used  to 
estimate  differences  in  effectiveness  between  iwo  types  of  MAD 
gear. 
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•  History  Reanalysis  where  one  run  was  made  unbiased  and  the  second 
generated  using  weighting  factors  on  the  histories  from  the  first 

to  correct  for  differences  in  MAD  gear.  This  was  also  used  for 
estimating  differences  in  effectiveness  between  two  types  of  MAD 
gear. 

•  Importance  Sampling  performed  with  an  importance  function 
weighted  to  generate  correlated  samples  for  the  two  types  of  MAD 
gear.  Again,  this  was  used  for  estimating  differences  in  effective¬ 
ness  between  two  types  of  MAD  gear. 

As  a  basis  for  quantitatively  characterizing  the  calculatior.al  efficiency 
of  variance  reduction  over  straightforward  or  crude  sampling  (i.  e. ,  with  no 
variance  reduction)  the  following  definition  was  used  throughout 

f  _  ,  variance  with  crude  sampling _ . 

"  Variance  with  the  variance  reduction  technique' 

,  APAIR  running  time  with  crude  sampling  . 

'  APAIR  running  time  with  variance  reduction' 

Crude  sampling  represents  the  procedure  currently  used  in  APAIR.  Using 
the  above  definition  for  efficiency,  a  value  of  e  =  2  for  application  of  a 
variance  reduction  technique  implies  a  reduction  in  computer  time  by  1/2 
to  achieve  the  same  variance  as  would  be  obtained  with  crude  sampling. 

The  rationale  for  using  this  as  an  efficiency  factor  is  discussed  in  Ref.  3. 

It  should  be  recognized  that  the  efficiency,  e ,  as  defined  above  is  a.  ran¬ 
dom  variable  since  the-  variances  with  and  without  variance  reduction  are  esti¬ 
mates  generated  from  the  statistics  of  the  simulation  and  not  theoretical  analyses. 
A  batching  method  was  used  in  estimating  most  of  these  variances.  The  proce¬ 
dure  is  presented  in  Ref.  3  and  will  not  be  detailed  here.  It  is  important  to 
recognize,  however,  that  the  efficiencies  presented  are  estimates  and,  there¬ 
fore,  are  subject  to  a  certain  variability.  In  fact,  the  variance  estimates  are 
second  order  quantities  and  thus  this  variability  is  generally  much  greater  than 
that  encountered  for  the  primary  quantity  whose  variance  is  being  estimated. 
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From  theoretical  arguments  the  error  in  the  efficiency  figure  could  easily  be 
as  large  as  —  45%.  In  some  of  the  cases  presented  in  this  report,  there  are 
large  correlations  between  the  crude  simulation  and  the  variance  reduced  tech¬ 
nique  and  thus  the  variability  of  the  efficiency  figure  is  much  less  than  this 
theoretical  maximum.  Since  efficiency  is  so  difficult  to  calculate  with  any 
accuracy,  calculations  were  made  for  five  similar  parameters— kill  probability, 
mean  time  to  kill,  number  of  stores  used  per  kill  (for  three  types  of  stores: 
torpedoes  and  two  kinds  of  sonobuoys) .  In  general  the  variance  reduction 
efficiencies  for  these  parameters  should  not  be  greatly  different,  therefore, 
the  spread  in  efficiency  values  for  these  five  parameters  should  give  a  rough 
idea  of  the  variability  of  the  efficiency  figure  for  each  variance  reduction 
technique.  As  an  improved  estimate  of  the  efficiency  of  each  technique,  a 
simple  average  of  the  five  efficiency  figures  is  presented  in  the  tables. 

The  initial  variance  reduction  technique  implemented  was  statistical 
estimation.  Implementation  involved  patching  into  the  report  generator 
portions  of  the  program  to  call  subroutines  designed  to  score  results  using 
statistical  estimation.  These  routines  contained  the  batching  needed  for 
estimation  of  variances,  print  statements  for  the  results,  etc.  To  minimize 
the  programming  changes  needed  to  implement  subsequent  techniques,  the 
statistical  estimation  subroutines  were  left  in  the  program  to  calculate 
and  print  results.  Thus  each  of  the  subsequent  comparisons  that  were  done 
actually  contrasted  a  variance  reduction  technique  plus  statistical  estimation 
to  a  case  with  only  statistical  estimation.  However,  with  the  exception  of 
the  expected  value  technique  as  discussed  in  Section  3.3,  it  was  felt  that 
the  efficiency  factors  obtained  were  essentially  the  same  as  if  the  variance 
reduction  technique  alone  had  been  compared  to  crude  Monte  Carlo. 

Each  of  the  variance  reduction  techniques  used  here  have  been  described 
in  detail  in  Ref.  3  and  will  not  be  presented  to  that  extent  in  this  report. 


13 


However,  in  the  following  discussions,  the  techniques  and  the  results  of  the 
application  of  the  techniques  will  be  described  in  sufficient  detail  to  provide 
an  appreciation  for  the  steps  involved.  Before  proceeding,  a  brief  descrip¬ 
tion  of  the  general  type  of  APAIR  problem  considered  in  the  study  will  be 
presented. 

3. 1  APAIR  PROBLEM  DESCRIPTION 

APAIR  can  be  applied  to  a  variety  of  ASW  problems  involving  one 
aircraft  in  pursuit  of  a  submarine.  The  problem  addressed  here  is  con¬ 
sidered  to  be  a  typical  application  of  APAIR  in  that  it  is  designed  to  esti¬ 
mate  the  effectiveness  of  the  aircraft  and  sensors  in  finding  and  killing  a 
submarine. 

A  total  simulation  of  APAIR  is  comprised  of  a  series  of  independent 
histories.  A  single  history  in  the  simulation  proceeds  as  follows: 

The  initial  aircraft  location,  direction  and  speed  are  chosen.  The 
submarine  location,  bearing  and  speed  are  also  chosen  according  to  spe¬ 
cified  random  characteristics.  The  aircraft  proceeds  to  execute  certain 
tactics  with  the  objective  of  detecting  the  submarine. 

Once  the  submarine  is  detected,  the  aircraft  enters  a  localization 
phase  where  certain  maneuvers  are  performed  and  scnobucys  are  dropped 
in  a  specified  pattern.  In  addition,  Magnetic  Anomaly  Detection  (MAD) 
gear  is  used  in  the  localization  process.  Two  typical  problems,  differing 
only  in  the  range  of  detection  of  the  MAD  gear,  were  utilized  in  this  study. 

If  localization  of  the  submarine  has  been  achieved,  the  aircraft 
proceeds  to  drop  a  torpedo  which  may  or  may  not  result  in  a  kill.  If  no 
kill  is  realized,  the  aircraft  then  proceeds  through  further  localization  to 
attempt  another  torpedo  drop,  again  following  specified  maneuvers.  The 
history  can  be  terminated  in  several  ways  which  include  exceeding  a  time 

limit,  achieving  a  submarine  kill,  or  exhausting  aircraft  stores  (torpedoes 
or  sonobuoys). 
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Once  the  history  is  terminated,  by  any  event  whatever,  a  new 
history  is  initiated.  When  the  specified  number  of  histories  are  completed, 
final  results  are  tabulated  and  printed  out.  Among  the  parameters  esti¬ 
mated  by  APAIR  the  following  were  used  in  this  study  to  illustrate  the  effi¬ 
ciency  of  variance  reduction: 

•  Probability  of  a  submarine  kill 

•  Average  time  used  to  achieve  a  submarine  kill 

•  Number  of  stores  used  per  submarine  kill  for  three  types  of 
stores.  One  type  was  torpedoes  and  the  other  two  were  two 
kinds  of  sonobuoys.  (As  the  exact  identification  of  the  sono- 
buoy  types  was  not  clear  from  our  documentation  and  irrele¬ 
vant  to  this  study,  they  are  merely  labelled  type  A  and  type  B 
in  this  report.  A  third  type  of  sonobuoy  was,  for  the  tactics 
in  our  sample  problem,  not  a  random  variable  and  thus  was 
not  included) 

Improvement  in  the  estimate  of  these  basic  parameters  was  the 
objective  in  using  the  variance  reduction  techniques  considered  here.  Of 
course,  such  an  improvement  can  also  be  achieved  by  increasing  the  num¬ 
ber  of  histories,  although  the  running  time  can  become  prohibitively  long. 
For  example,  the  above  problem  required  approximately  15  minutes  on  a 
Univac  1108  to  generate  100  histories.  The  desirability  for  achieving 
variance  reduction  is  therefore  obvious. 
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3.  2  APPLICATION  OF  STATISTICAL  ESTIMATION 

In  the  course  of  an  individual  history,  the  aircraft  will  drop  a  torpedo 
when  it  thinks  it  has  determined  the  submarine 's  location  and  heading.  The 
actual  submarine  speed,  aspect,  and  range  are  used  to  determine,  from  input 
tables,  the  probability,  p^,  that  the  torpedo  will  destroy  the  submarine.  A 
random  number,  Ry,  from  a  uniform  distribution  is  generated  and  compared 
to  p^.  If  R  <  Pr,,  a  kill  is  scored  and  the  history  is  terminated.  If  R  >  p„, 
no  kill  occurs  and  no  scoring  is  done;  the  game  continues  with  more  localization 
and  possibly,  another  attack.  If  other  torpedoes  are  dropped,  a  similar  proce¬ 
dure  is  used  to  determine  if  a  kill  is  scored  later  in  the  game. 

At  the  completion  of  the  simulation,  the  probability  of  submarine  kill  is 
estimated  as 

PK  =  n/N 

where 

n  =  number  of  kills  scored  and 
N  =  number  of  histories  run. 

For  estimating  the  time  to  kill,  APAIR  currently  uses 

%  =  x/n  2  TKi 
i-1 

where  n  is  the  number  of  kills  scored  and  TKi  ;  i  =  l, . . . ,  n  the  time  taken 

th  ^ 

to  effect  a  kill  in  the  i  history  which  ended  in  a  kill.  Similar  procedures  are 
used  for  the  other  parameters  (number  of  torpedoes  and  sonobuoys  used  per 
kill)  estimated. 
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In  the  statistical  estimation  technique  for  variance  reduction,  the  scoring 
was  changed  so  that  actual  kills  are  not  scored,  but  rather  the  expected  value 
of  kill,  Pr,,  was  scored  for  all  torpedo  drops  regardless  of  whether  or  not  a  kill 

IV 

was  achieved.  However,  the  simulation  game  itself  was  not  modified  in  any 
way.  That  is,  a  random  number,  was  still  used  at  each  torpedo  drop 
to  determine  whether  the  history  was  terminated  by  a  kill  or  more  maneuvers 
took  place.  The  outcome  of  this  random  choice  did  not  affect  the  scoring  of  p^.. 


Specifically,  define 


probability  of  kill  for  the 


.th 

3 


torpedo  dropped  during  history 


i. 


n.  =  number  of  times  a  torpedo  was  dropped  during  history  i . 
Then,  the  total  score  generated  for  history  i  if. 


i  =  l> 


,N 


and  the  estimate  for  the  probability  of  kill  for  the  entire  simulation  (N  histories) 
is 

N  ni 

PK  =  ly/N  Xv  2  pKij 

i=l  i=l 
* 

In  the  case  of  the  remaining  parameters,  a  similar  approach  was  taken. 
For  example  if 

tTi 

T™.  =  time  to  torpedo  detonation  for  the  j  torpedo  dropped  in 
3  history  i 

then  the  estimate  for  time  to  kill  is  given  by 
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N  “i 

V  l/N  Z  Z  PKijTKii 

1=1  ]=1 


Similar  egressions  were  used  for  the  remainder  of  the  parameters  being 
estimated. 

The  computational  time  per  history  is  very  slightly  increased  using 
this  technique  since  the  same  game  is  still  being  played  but  there  is  a  little 
more  bookkeeping  in  the  scoring.  However,  this  is  offset  by  the  resulting 
variance  reduction. 

In  demonstrating  the  statistical  estimation  technique,  both  the  crude 
Monte  Carlo  estimate  (counting  actual  kills)  and  the  statistical  estimate  (sum¬ 
ming  over  pK  values)  were  calculated  in  the  same  run  and,  therefore,  used 
the  same  histories.  This  produces  a  high  degree  of  correlation  between  the 
crude  and  the  statistical  estimation  results  which  reduces  the  variance  of  the 
efficiency  figure.  Two  problems  were  run  using  MAD  gear  having  long  and 
short  range  detection  capabilities  respectively.  The  variances  obtained  with 
statistical  estimation  and  with  crude  Monte  Carlo  are  show  in  Tables  3. 1  and 
3. 2  along  with  the  resulting  efficiency  factor  for  the  use  of  this  variance  reduc¬ 
tion  technique.  The  sample  variances  were  estimated  using  the  statistical 
techniques  described  in  Ref.  3  to  obtain  the  variance  results  indicated.  The 
actual  estimated  values  of  the  parameters  are  not  presented  since  they  are  not 
considered  to  be  germane  to  comparison  of  the  efficiencies.  Therefore,  only 
the  variances  in  the  two  cases  are  shown  along  with  the  efficiencies  obtained 
in  estimating  the  parameters  with  variance  reduction. 

The  efficiencies  obtained  varied  from  1. 00  (implying  no  improvement) 
to  as  high  as  1. 50  (implying  a  factor  of  1, 5  reduction  in  running  time).  However, 
these  extreme  values  probably  represent  statistical  fluctuations  in  the  variance 
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TABLE  3. 1 

stimation  Technique  For  the  Short  Range  MAD  Gear 
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TABLE  3.  2 

Variance  Reduction  Using  the  Statistical  Estimation  Technique  For  the  Long  Range  MAD  Gear 


Average  Efficiency  Improvement  s  1.  23 

Ratio  of  time  increase  with  variance  reduction  s  1.01 

Statistical  Estimation:  Rather  than  score  actual  kills,  an  expected  value  of  kill  is  scored  at  each 
torpedo  drop  regardless  of  v/hether  or  not  the  torpedo  scored  a  kill.  The  game  is  not  modified. 


estimates  rather  than  real  efficiency  improvement.  Average  efficiencies  of 
1. 15  for  the  short  MAD  problem  and  1. 23  for  the  long  were  achieved.  Thus, 
it  can  be  seen  that  use  of  statistical  estimation  could  have  made  roughly  a  19% 
improvement  in  APAIR  run  times.  Use  of  statistical  estimation  is  justified  since  it 
involves  a  very  trivial  modification  to  the  program. 

3. 3  EXPECTED  VALUE  TECHNIQUE  FOR  ESTIMATING  EFFECTIVENESS 

OF  MULTIPLE  TORPEDO  DROPS 

The  expected  value  technique  differs  from  statistical  estimation  in  that 
the  actual  game  or  simulation  being  played,  as  opposed  to  just  the  scoring  tech¬ 
nique,  is  changed  to  replace  a  random  choice  with  an  expected  value  for  the  out¬ 
come  of  that  choice.  For  example,  instead  of  generating  Ry  to  test  against 
PK  with  the  choice  of  either  killing  or  missing  the  submarine,  the  submarine  is 
given  an  initial  "weight"  of  1. 0.  If  the  first  torpedo  dropped  has  a  pK  of  .  80, 
then  80°?*  of  the  submarine  is  deemed  killed  but  20%  survives  and  the  weight  is 
reduced  to  .2.  If  a  second  torpedo  is  dropped  which  has  a  pK  of  .  50,  then  half 
of  the  remaining  weight,  or  .  1,  is  killed,  while  a  weight  of  .  1  continues  to 
survive.  The  history  is  never  ended  due  to  a  submarine  kill  but  continues  until 
some  other  limit,  such  as  using  up  the  mission  time  or  using  up  the  stores  of 
torpedoes  or  sonobuoys,  stops  the  history. 

Such  a  technique  could  be  useful  in  studies  such  as  determining  the 
effectiveness  of  the  number  of  torpedoes  carried  on  an  aircraft  or  the  worth¬ 
iness  of  the  tactics  for  relocalization  and  reattack  following  a  miss  by  the 
first  torpedo.  In  these  cas^s  one  is  not  so  much  interested  in  the  overall  kill 
probability  as  in  the  kills  scored  by  the  second,  third,  etc.  torpedoes.  In  a 
crude  Monte  Carlo  most  of  the  kills  are  made  by  the  first  torpedo  and  the  his¬ 
tory  ends  there.  These  histories  add  nothing  to  the  knowledge  of  tactics  in¬ 
volving  reattack  or  to  the  kill  value  of  the  second  torpedo,  but  simply  con¬ 
stitute  wasted  time  towards  calculating  the  items  of  importance.  In  fact,  what 
is  worse,  these  histories  add  variance  to  the  overall  kill  probability.  It  would 
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be  prohibitively  expensive  to  determine  the  value  of  reattacks  from  a  crude 
Monte  Carlo  due  to  the  large  proportion  of  the  running  time  being  spent  on 
histories  of  no  value  to  this  parameter. 

Using  the  expected  value  technique  outlined  above,  most  histories  will 
contribute  to  knowledge  concerning  the  second  and  third  torpedoes  as  the  history 
will  always  continue  after  the  first  torpedo  drop.  For  example,  one  of  the  crude 
Monte  Carlo  problems  run  had  (out  of  100  histories)  78  cases  in  which  the  first 
torpedo  was  dropped,  10  in  which  the  second  was  dropped,  and  only  one  history 
where  the  third  torpedo  was  dropped.  Using  the  expected  value  technique,  the 
same  problem  had  79  histories  with  one  torpedo  drop,  65  histories  in  which  the 
second  torpedo  was  dropped,  and  43  histories  of  a  third  torpedo  drop. 


Estimation  in  the  expected  value  case  was  done  as  follows. 

tii  tii 

probability  for  the  j  n  torpedo  in  the  i  history  and  w. . 
weight  at  the  time  that  torpedo  was  dropped,  then 


If  pjQj  is  the  kill 
is  the  submarine 


N 

PKj  =  1/N  PKijwij 
i=l 

is  the  kill  probability  of  the  torpedo.  Likewise,  the  time  to  kill  by  the 
torpedo  is  given  by 


N 

1  Kj  =  /f~  Zj  TKijpKijwij 
PKj  i=i 

where  TR. .  is  the  time  of  detonation  of  the  j  torpedo  on  the  i  history. 

A  similar  formula  is  used  to  estimate  the  numbers  of  sonobuoys  dropped  per 
th 

kill  by  the  j  torpedo.  (Obviously  the  number  of  torpedoes  used  is  constant 
in  this  case,  so  this  parameter  was  omitted. ) 
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For  this  technique,  two  sample  problems  were  run;  one  with  the  long- 
range  MAD  gear  and  one  with  the  short-range  MAD  gear.  In  each  case  APAIR 
was  run  twice,  with  and  without  the  expected  value  technique.  To  reduce  the 
variance  of  the  efficiency  factors,  the  histories  in  the  two  runs  were  correlated 
using  the  technique  described  in  Section  3.  6.  This  kept  the  histories  in  each 
run  identical  through  the  first  torpedo  drop.  The  resulting  efficiencies  for  the 
second  and  third  torpedo  drops  for  botli  cases  are  shown  in  Table  3.  3.  There 
were  no  examples  of  a  fourth  torpedo  drop  in  any  of  the  runs.  In  the  short  MAD 
case,  there  were  no  examples  of  a  third  torpedo  drop  in  the  crude  Monte  Carlo 
run,  so  the  efficiency  of  the  variance  reduction  technique  is  theoretically  in¬ 
finite  in  this  case.  However,  using  the  pK  estimate  from  the  expected  value 
run,  it  was  possible  to  estimate  the  number  of  crude  Monte  Carlo  histories  that 
would  be  necessary  to  get  similar  statistics;  this  led  to  the  efficiency  factor  of 

a 

10  for  pK  shown  for  the  third  torpedo  in  the  short  MAD  case. 

The  running  times  for  the  crude  and  expected  value  calculations  are 
show  in  Table  3. 4.  As  anticipated,  the  expected  value  histories  took  much 
longer  to  run  because  they  did  not  stop  at  the  first  torpedo  drop  (as  most  of 
the  crude  Monte  Carlo  histories  did)  but  went  on  to  simulate  a  second  and  a 
third  torpedo  drop.  This  extra  running  time  is  used  to  generate  information 
concerning  the  parameters  of  interest,  i.  e. ,  the  kills  made  by  second  and  third 
torpedo  drops,  and,  therefore,  the  overall  efficiency  is  much  improved  as 
Table  3. 3  shows.  It  is  instructive  to  consider  the  efficiency  for  the  first 
torpedo  drop.  As  the  histories  in  the  two  runs  were  identical  through  the  first 
torpedo  drop,  the  variances  are  identical  for  the  first  torpedo  parameters. 

Due  to  the  increased  running  time,  the  first  torpedo  efficiencies  will  be  lower 
by  .  49  (for  the  short  MAD)  and  .  36  (for  the  long  MAD).  The  extra  running  time 
used  in  calculating  second  and  third  torpedo  drops  is  wasted  as  far  as  first 
torpedo  drops  is  considered  and  this  lowers  the  efficiency.  This  illustrates 
a  common  point  of  variance  reduction:  any  technique  which  reduces  variance 
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TABLE  3. 4 

Running  Times  for  Crude  Monte  Carlo  and  Expected  Value  Calculations 


Problem 

Running  Time  for 
'Crude*  Monte  Carlo 

Running  Time  for 
Expected  Value  Technique 

Ratio 

of 

Times 

Problem  1: 

Long  MAD 

664  sec 

1866  sec 

2.81 

Problem  2: 

Short  MAD 

822  sec 

1666  sec 

2.03 
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for  one  parameter  will  increase  variance  for  some  other  parameter.  Variance 
reduction  techniques  must  be  carefully  tailored  to  the  parameters  of  importance, 
in  this  case  the  kills  involving  two  and  three  torpedo  drops. 

For  ease  in  making  programming  changes  to  calculate  the  parameters 
separately  by  torpedo  drop,  base  runs  with  the  current  APAIR,  representing 
crude  Monte  Carlo,  were  not  made.  Runs  using  the  expected  value  technique 
then  provided  efficiency  factors  for  an  expected  value/statistical  estimation 
comparison.  The  average  efficiency  of  the  statistical  estimation/crude  Monte 
Carlo  comparison  in  Section  3. 2  was  determined  to  be  a  factor  of  1. 2.  Multi¬ 
plying  the  efficiencies  from  the  expected  value/statistical  estimation  compari¬ 
sons  by  1. 2  generated  the  figures  presented  in  Table  3. 3  as  the  efficiency  of 
expected  value  versus  crude  Monte  Carlo.  This  combination  of  efficiency 
factors  is  justified  because  the  expected  value  technique  necessarily  incorpo¬ 
rates  statistical  estimation  scoring;  once  the  kill/miss  decision  has  been  re¬ 
moved  from  the  game,  the  scoring  must  be  by  the  p  *s  for  each  torpedo  drop. 
Thus  there  is  no  possible  "expected  value  without  statistical  estimation"  com¬ 
parison  to  crude  Monfe  Carlo  but  only  the  combined  efficiency. 

3. 4  SYSTEMATIC  SAMPLING  OF  INITIAL  SUBMARINE  POSITION 

Systematic  sampling  is  a  variance  reduction  technique  that  usually  finds 
application  in  selecting  initial  or  starting  values  for  a  random  variable.  Poten¬ 
tial  applications  in  APAIR  include  initial  submarine  or  aircraft  bearing  and 
location.  Sampling  in  a  systematic  manner  essentially  serves  to  reduce  the 
contribution  to  the  variance  coming  from  the  random  variables  being  systema¬ 
tically  sampled. 

To  demonstrate  this  technique  in  APAIR,  problems  were  run  with  the 
aircraft  initial  location  and  both  the  aircraft  and  submarine  initial  bearings 
fixed.  This  left  the  submarine  starting  location  as  the  variable  which  was 
systematically  sampled  as  shown  in  Fig.  3. 1.  The  aircraft  was  initially 
located  at  the  origin  while  the  submarine  starting  position  was  uniformly 
distributed  along  the  y-axis  (north)  between  0  and  L.  Both  the  aircraft  and 
the  submarine  were  initially  moving  east  as  shown. 
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N 


Location 

Fig.  3. 1.  Starting  Positions  for  Systematic  Sampling 
Demonstration. 


Systematic  sampling  was  applied  in  two  different  ways.  The  first 
was  used  on  the  long  MAD  gear  problem  and  was  implemented  as  follows: 

To  obtain  starting  positions  for  the  first  10  histories,  a  random  num¬ 
ber  Ry  was  selected  from  U(0, 1)  and  ten  initial  submarine  positions  were 
located  at 

Vj  =  LRU/10 

y„  =  L/10  +  LR../10 

a  u 

y3  =  2L/10  +  LRU/10 


y1Q  =  9L/10  +  LRu/l0 

The  histories  run  using  these  initial  starting  conditions  constituted  the  first 
batch.  Then  another  series  of  ten  y's  was  generated  by  selecting  a  second 
random  number  from  U(0, 1),  and  the  second  batch  was  run.  The  process 
was  continued  until  a  total  of  10  batches  of  results  (or  100  histories)  were 
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obtained.  The  results  of  these  simulations  were  used  to  estimate  the  parame¬ 
ters  and  sample  variances  using  batched  estimators.  That  is ,  an  estimate 
for  Pj^  was  obtained  for  each  batch  by  the  usual  methods.  These  batch  esti¬ 
mates  may  be  denoted  Pjq,  *  *  *  ,  A  final  estimate  was  obtained  from 


K10 


10 


■K 


■  1/10 E 


•Ki 


i=l 


and  the  sample  variance  was  estimated  from 

x2 


1C 

E 

i=l 


S2  =  1/9  £  (pK1  -  5k): 


The  results  are  shown  in  Table  3. 5  which  summarizes  the  variances  obtained 
with  and  without  systematic  sampling.  It  can  be  seen  that  the  results  are 
rather  mixed,  and  in  some  cases  a  worse  result  was  obtained  using  systematic 
sampling.  Most  of  this  variation  is  strictly  statistical  fluctuation  due  to  the 
unavoidably  large  variances  in  the  efficiency  estimates.  That  this  is  the  case 
may  be  seen  from  the  efficiencies  which  are  less  than  1. 0.  Theoretically, 
systematic  sampling  should  always  reduce  variance  and  efficiencies  should 
always  be  1.0  or  greater.  However,  if  efficiencies  are  close  to  1,  it  is  easy 
to  get  estimates  which  are  just  below  1.  0. 

Some  of  the  variation  in  the  systematic  sampling  efficiencies  may  also 
represent  a  variation  in  how  sensitive  a  parameter  is  to  the  submarine  starting 
location.  The  probability  of  kill  (eventually,  after  enough  localization)  should 
not  be  as  dependent  on  the  submarine’s  starting  distance  as  the  time  taken  to 
localize  and  kill  or  the  number  of  sonobuoys  used  in  localization.  Any  parame¬ 
ter  which  is  not  sensitive  to  submarine  starting  position  will  have  a  variance 
which  is  likewise  not  sensitive  to  reduction  of  variance  in  selecting  the  sub¬ 
marine  starting  position. 
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Variance  Reduction  Using  Systematic  Sampling 


The  average  efficiency  gain  for  systematic  sampling  in  the  long 
MAD  problem  was  a  factor  of  1. 34. 

A  second  type  of  systematic  sampling  was  used  in  the  problem 
with  the  short  MAD  gear.  In  this  case,  a  new  random  number  was  used 
each  time  a  new  starting  position  for  the  submarine  was  selected.  Ten 
strata  were  used  with  the  starting  position  limited  to  fall  between  ranges 
and  Lg.  That  is,  the  first  ten  starting  positions  were  selected  using 


(Lo  -  LJ 

yl  "  L1  +  - TO -  Rul 


<L2  -  Ll> 


y2  ~  L1  +  - H5 -  (IS  + 


y*  = 


(L9  -  L  ) 

L<  +  (Ru„  +  2.) 


3  ~  n  T  — nr 


'10 


L1  +  -T-  <So  +  9,1 


where  RUl,  *  •  • ,  Ruin  are  random  samples  from  U(0, 1). 


The  samples  in  this  case  were  batched  as  before  and  the  results 
presented  in  Table  3. 6  were  obtained.  These  results  are  seen  to  be  quite 
similar  to  those  obtained  for  the  long  MAD  gear.  Again,  a  large  variance 
in  the  efficiency  is  apparent.  Also,  it  appears  that  several  of  the  param¬ 
eters  were  insensitive  to  initial  submarine  starting  position.  The  average 
efficiency  is  1. 28  with  an  overall  uncertainty  of  about  20%  expected  in  the 
efficiency  estimate. 
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TABLE  3.6 

Variance  Reduction  Using  Systematic  Sampling  For  The  Short  Range  MAD  Gear 


3. 5  ANTITHETIC  VARIATES  FOR  SAMPLING  INITIAL  SUBMARINE 

POSITION 

The  final  variance  reduction  technique  applied  to  selection  of  the 
submarine  starting  position  was  the  use  of  antithetic  variates.  This  was 
performed  for  the  same  two  situations  described  for  systematic  sampling 
in  the  previous  sections. 

In  its  simplest  application,  use  of  antithetic  variates  seeks  to  generate 
negatively  correlated  samples  by  selecting  two  values  x* ,  x"  of  the  random 
variable  from  the  distribution  f(x)  using 


f(x)dx 


and 


!-Ru  *  J. 


f(x)dx 


where  Ry  is  a  random  number  selected  from  U\\  1).  The  values  of 
x'  and  x"  are  clearly  correlated  since  they  have  been  generated  by  the 
same  random  number  R^.  Also,  x*  and  x"  are  negatively  correlated 
since  when  x*  is  large,  z"  will  be  smalL 


In  the  application  here,  pairs  of  initial  starting  positions  for  the  sub¬ 
marine  were  selected  according  to  the  above  formula!*  Mi.  Thus,  when  one 
submarine  starting  position  was  selected  far  from  the  aircraft,  a  position 
close  to  the  aircraft  was  also  selected  fc  -he  nex*  history. 


In  the  first  problem,  (i.  e. ,  where  the  long  range  MAD  gear  was  used 
and  the  submarine  was  located  between  0  and  L),  the  pairs  of  starting 
positions  were  obtained  using 


y..  =  R  L 
u 

y2  *  -  Hu)L 
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and  two  histories  were  run. 


In  the  second  problem  (i.  e. ,  with  the  short  range  MAD  gear  and 
where  the  submarine  was  located  between  Lj  and  L2)  the  pairs  of  starting 
positions  were  obtained  using 

>'l  =  Lj  +  (I<2  -  1.^ 

H  -  L2  +  <Li  -  L2)Ru 

and  the  two  histories  corresponding  to  these  initial  starting  positions  were 
run. 

In  both  cases,  batching  was  performed  to  obtain  estimates  for  the 
variances.  The  results  of  the  analyses  are  summarized  in  Tables  3. 7  and 
3. 8  respectively. 

As  was  the  case  with  systematic  sampling,  there  was  a  wide  variation 
in  efficiency  and  one  variable  gave  worse  results  with  the  antithetic  variates 
than  with  crude  sampling,  indicating  as  before,  a  large  variance  for  the 
efficiency  estimates  and,  possibly,  several  parameters  which  were  not 
sensitive  to  the  initial  starting  yalues.  Although  it  is  theoretically  possible 
that  the  use  of  antithetic  variates  could  give  worse  results  than  crude 
sampling,  it  was  not  expected  here  and  the  single  value  less  than  X.  0  is  prob¬ 
ably  a  low  estimate  for  an  efficiency  just  above  X.  0.  In  any  event,  average 
efficiencies  of  1. 12  and  1. 37  were  obtained  in  the  two  problems.  An  ‘‘rror 
of  about  ±10%  is  expected  in  these  efficiencies. 

3.  6  CORRELATED  SAMPLING  FOR  ESTIMATING  DIFFERENCES  IN  MAD 

GEAR  EFFECTIVENESS 

Correlated  sampling  is  a  procedure  that  can  be  used  to  reduce  variance 
in  Monte  Carlo  simulation  in  the  following  general  situations: 

•  The  effect  of  a  perturbation  to  a  known  problem  is  to  be 
determined. 
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TABLE  3.  7 

Variance  Reduction  Using  Antithetic  Variates  For  The  Long  Range  MAD  Gear 


1 


TABLE  3.  8 

Variance  Reduction  Using  Antithetic  Variates  For  The  Short  Range  MAD  Gear 


•  The  difference  between  estimated  parameters  in  two  problems 
having  similar  characteristics  is  to  be  calculated. 

•  A  parametric  study  of  several  similar  problems  is  to  be 
performed. 

Such  situations  occur  very  frequently.  In  fact,  in  most  studies  the  most  im¬ 
portant  result  to  be  investigated  is  the  change  in  response  of  the  system  as  a 
problem  characteristic  is  varied.  As  will  be  seen,  the  payoff  from  various 
types  of  correlation  in  sampling  can  be  very  high. 

Use  of  the  APAIR  model  could  easily  involve  problems  having  one  or 
more  of  these  characteristics.  For  example,  the  sensitivity  to  a  range  of 
tactics  presents  a  potential  situation  where  correlated  sampling  could  provide 
substantial  improvements  in  efficiency. 

The  problem  selected  for  demonstration  of  correlated  sampling  in¬ 
volved  the  short  range  and  long  range  MAD  detector  cases  discussed  pre¬ 
viously.  The  main  parameter  of  interest  to  be  calculated  was  the  difference 
between  these  two  cases  in  the  parameters  used  in  prior  sections.  That  is, 
differences  in: 

•  Probability  of  submarine  kill 

•  Time  to  submarine  kill 

•  Number  of  torpedos  used  per  kill 

•  Number  of  sonobuoys  of  Types  A  and  B  used  per  kill. 

The  only  difference  in  problem  characteristics  in  the  two  cases  was 
a  difference  in  MAD  detection  capabilities.  This  was  expressed  in  a  function 
relating  probability  of  detection  to  range  of  target  from  the  aircraft,  as  shown 
in  Fig.  3. 2. 
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To  demonstrate  the  concept  of  correlated  sampling  for  this  applica¬ 
tion,  we  define 

Ps  =  probability  of  submarine  kill  using  the  short  range  MAD  detectors 

p^  =  probability  of  submarine  kill  using  the  long  range  MAD  detectors 

The  problem  of  interest  is  to  perform  Monte  Carlo  simulations  to  estimate 
the  difference 

A  ■  PL-ps 

in  the  case  of  probability  of  submarine  kill.  Similar  definitions  apply  to  the 
remainder  of  the  parameters  of  interest. 

The  crude  Monte  Carlo  approach  would  first  estimate  p„  (denoted  by 

A  A  b 

Pg  )  and  then  pL  (denoted  by  p^)  using  another  set  of  independent  histories. 
Then  the  difference  A  is  estimated  using 

A  A  A 

A'  =  ?i-V's 
The  variance  in  is 

o2(A')  ff2(Pg)  +  cr2(p^) 

for  the  case  of  independent  histories  in  the  Pg  and  pL  estimations.  Suppose, 
however,  that  nositivelv  correlated  estimates  for  o~  (sav  d_)  and.  for  t>, 

(say  $L)  were  used  to  estimate  A.  Then  for 

A  A  A 

A  =  PL'PS 

A 

the  variance  in  A  will  be  given  by 

cr"(A)  =  ff2(PL)  +  cr2(pg)  -  2  cov  (pg,PL) 
where  cov  (Pg,PL)  is  the  covariance  between  pg  and  pL*  Since  pg  and 
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pT  are  positively  correlated,  then 

Ju 

cov  (pg,p^)  2:  0,  and  hence 

or^(A)  ^  cr^(A')  . 

Thus,  the  objective  of  correlated  sampling  is  to  develop  a  sampling 

A  A 

strategy  that  will  correlate  the  estimators  pg  and  p^.  This  technique  can 
be  extremely  powerful  when  small  perturbations  of  problems  are  to  be  studied 
since  the  correlations  induced  will  tend  to  emphasize  differences  in  the  prob¬ 
lems  due  to  the  perturbation  rather  than  differences  due  to  statistical  fluctua¬ 
tions  which  are  usually  the  controlling  differences  in  cases  with  independent 
histories. 

The  above  arguments  for  improvement  in  the  variance  using  correlated 
sampling  for  the  probability  of  kill  would  also  apply  directly  to  estimating  the 
differences  in  other  ASW  parameters. 

Correlation  can  be  accomplished  in  several  ways.  For  example,  the 
short  and  long  range  MAD  problems  could  be  simulated  independently  except 
the  same  random  number  could  be  used  in  determining  the  outcome  once  the 
probability  of  detection  had  been  obtained  from  the  curves  in  Fig.  3. 2.  Thus, 
correlation  between  the  two  results  would  exist.  Another  way,  and  the  one 
that  was  used  here,  is  to  control  the  random  numbers  in  the  two  simulations, 
by  using  the  same  random  numbers  in  the  two  problems  until  a  diffemnce  in 
detection  occurs  in  the  problems  due  to  the  difference  in  the  MAD  gear.  Two 
separate  runs  were  made,  but  corresponding  histories  in  the  two  simulations 
were  made  to  start  at  the  same  point  in  the  random  number  sequence.  Thus,  the 
histories  would  be  identical  up  to  the  point  where  the  difference  in  MAD  gear 
resulted  in  different  decisions.  At  the  time  the  detectors  came  into  play,  the  de 
tecticn  outcome  was  selected  in  each  case  from  the  same  random  number.  Subse 
quently,  the  histories  continued  independently  until  the  end  of  the  game.  More 
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correlation  could  have  been  introduced  by  subsequently  using  identical  random 
numbers  wherever  the  problem  logic  allowed. 


The  program  changes  made  to  induce  this  much  correlation  were 
fairly  simple.  Two  separate  random  number  generators  were  used  in  each 
simulation.  The  first  was  used  once  each  history  to  give  a  random  starting 
point  in  the  sequence  of  the  second  generator;  it  produced  the  same  sequence 
of  starting  points  in  both  simulations.  The  second  generator  was  used  in  the 
history  to  obtain  random  numbers  for  the  simulation  process.  By  starting  at  the 
same  point  in  both  cases,  identical  histories  will  be  generated  until  there  is  a 
difference  in  decisions  made  concerning  a  MAD  detection.  Even  though  the  first 
history  in  one  problem  might  use  more  random  numbers  than  in  the  other 
problem,  the  random  number  sequences  would  be  returned  to  the  same  point 
at  the  start  of  the  second  history  in  each  problem. 


To  calculate  the  variance  of  the  difference.  A,  in  the  two  cases,  it 

a 

was  necessary  to  obtain  ’batch’  values,  A  ,  for  the  difference  in  the  batch 
values  pLn  and  pgn  that  have  been  described  in  earlier  sections.  The  batch 
values  of  pgn,  p^n,  and  the  otner  parameters  were  written  on  temporary 
files  as  the  simulations  progressed.  Then  a  separate  small  program  was 
written  to  combine  these  files  and  calculate  the  batch  differences  in  the  vari¬ 
ous  parameters.  The  final  estimated  difference  was  the  average  of  the  batch 
differences  and  the  variance  of  the  difference  was  determined  from  the  spread 
of  the  sample  batch  differences.  Specifically,  by  grouping  the  100  histories 
of  the  long  MAD  simulation  into  batches  of  10,  one  calculates 

10 

PLi  =  1/10  PKi 
i=l 


20 


PL2  =  1/10  L  % 


i=ll 
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100 

PL10  =  1/10  S  PKi 

i=  91 

where  pLr  is  the  estimate  of  p^-  from  batch  n  and  P-^-  is  the  estimate  of 
PK  from  history  i.  Similar  formulas  were  used  to  obtain  batch  estimates  for 
the  other  parameters  in  the  long  and  short  simulations.  Then  the  batch  differ¬ 
ences  were  calculated; 

h  =  hi-  hi 


A  A  A 

A10  =  PL10“PS10  ' 

The  final  estimator  for  the  difference  in  probability  of  kill  is 
10 

4  =  1/10  ]T  4n 

n=l 

and  the  estimated  variance  is 


-2,;.  1 
a  M  =B 

10 

4Z 

m 

(A  -  A)2 
v  n 

— <105 

III 

10 

/  10  \2 
(1/10E  hi 

n=l 

n=l 

\  n=l  / 

Rather  than  make  two  additional  uncorrelated  runs  to  get  comparison  variance 
estimates  for  the  uncorrelated  or  crude  Monte  Carlo,  the  uncorrelated 
equation 

s2(4)  =  »2(ps)  + 


was  used  with 


aO  a 

(Pq) 


and 


The  results  of  the  calculations  are  shown  in  Table  3. 9.  The  greatest 
increase  in  efficiency  (a  factor  of  3)  was  found  for  the  probability  of  kill.  The 
average  improvement  over  all  parameters  was  almost  a.  factor  of  2.  Thus, 
the  running  time  for  the  same  variance  could  be  reduced  by  about  a  factor  of 
2  by  using  correlated  sampling. 

There  was  a  slight  increase  in  running  time  (although  the  vari¬ 
ance  was  substantially  reduced)  with  variance  reduction  since  some  additional 
bookkeeping  was  used  in  the  program.  The  running  time  could  have  been  re¬ 
duced  with  some  additional  effort  by  not  recalculating  v/ith  identical  random 
numbers  the  part  of  each  history  up  to  the  point  where  the  detection  came  into 
play,  but  by  simply  saving  the  results  of  one  case  for  application  to  the  other. 
However,  this  would  have  involved  more  extensive  computer  program  modifica¬ 
tions  than  were  warranted  here. 

Analysis  of  Correlated  Histories 

One  of  the  great  benefits  of  correlation  is  the  potential  it  provides  to 
gain  insight  into  understanding  simulation  problems.  For  example,  if  two 
highly  correlated  histories,  one  with  and  one  without  a  problem  perturbation 
were  available,  then,  most  of  the  time,  differences  observed  in  the  histories  will 
be  due  to  variations  in  the  problem  perturbation  rather  than  statistical  variations. 

The  problem  of  the  two  MAD  detectors  was  ideally  suited  for  demon¬ 
strating  the  possibility  of  analysis  of  correlated  histories  since  the  MAD  detectors 
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TABLE  3. 9 


Variance  Reduction  Using  Correlated  Sampling  For  Estimating  Differences 
Between  Short  and  Long  Range  MAD  Gear  Parameters 


Estimated  Difference 
in  Expected  Value 
(Long  Range  MAD  - 
Short  Range  MAD) 

Variance  With 
Crude  Sampling 

Variance  With 
Correlated  Sampling 

Efficiency 

Probability  of 
Submarine  Kill 

2.96  x  10 "3 

9.6  x  10“4 

3.0 

Time  to 

Submarine  Kill 

15.0 

9.2 

1.58 

Number  of  Torpedos 
per  Kill 

2. 2  x  10 "3 

1.6  x  10  "3 

1.33 

Sonobuoys  Used  per 
Kill  -  Type  A 

0.166 

0.068 

2.40 

Sonobuoys  Used  per 
Kill  -  Type  B 

0. 174 

0.144 

1.17 

Average  Efficiency  Improvement  s- 1. 9 

Ratio  of  time  increase  with  variance  reduction  s:  1. 03 

Correlated  Sampling:  Two  separate  runs  were  made.  Histories  were  identical 
up  to  point  where  difference  in  MAD  gear  lead  to  a  divergence  in  histories. 


1.0 
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were  sufficiently  similar  so  that  variations  on  estimates  for  the  kill 
probability  could  be  lost  in  the  statistics  of  the  problem.  More  specifically, 
consider  the  sample  series  of  single  histories  shown  in  Table  3. 10.  These 
histories  were  obtained  from  runs  correlated  in  the  manner  described  above. 

There  are  three  types  of  situations  that  can  be  identified  in  Table 
3. 10.  First,  there  are  a  large  number  of  histories  where  the  short  and  long 
range  gear  give  the  same  result.  This  situation  should  be  expected  since  the 
gear  is  similar  and  the  histories  are  highly  correlated.  That  is,  when  a  short 
range  detection  occurs  with  the  short  MAD,  it  will  also  occur  with  the  long  MAD. 
Also  when  no  detection  occurs  for  the  long  MAD  at  long  ranges,  none  will  occur 
for  the  short  MAD. 

At  intermediate  ranges  where  the  two  MAD  curves  differ,  this,  of 
course,  is  not  true.  This  is  manifested  by  the  second  situation  (history  8) 
where  the  long  range  MAD  detected  the  submarine  and  effected  a  kill  and  the 
short  range  MAD  didn't,  with  no  kill  as  a  result.  This  type  of  history  was  also 
expected. 

Of  most  interest  is  the  third  situation  where,  in  histories  i  and  9,  the 
long  range  MAD  gear  produces  a  lower  kill  probability  than  the  short  MAD 
gear  although  both  types  detected  the  submarine.  This  result  was  quite  un¬ 
expected  and  would  tend  to  indicate  there  are  considerations  other  than  varia¬ 
tions  in  MAD  gear  which  make  the  two  problems  different.  For  example,  the 
tactics  used  might  be  appropriate  to  the  short  MAD  detector  but  might  be 
'trigger  happy  *  when  used  with  the  long  MAD  gear  resulting  in  premature  torpedo 
drops  with  a  lower  kill  probability.  Had  the  histories  not  been  correlated, 
it  would  have  been  difficult,  due  to  the  statistical  variation  from  history  to 
history,  to  notice  that  such  events  were  occurring.  In  this  manner,  therefore, 
correlation  can  identify  which  histories  should  be  examined  in  detail  to  provide 
more  insight  into  the  problem. 
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Another  application  of  correlated  sampling  is  presented  in  the  follow¬ 
ing  section. 

3. 7  HISTORY  REANALYSIS  FOR  ESTIMATING  DIFFERENCES  IN  MAD 

GEAR  EFFECTIVENESS 

In  the  second  application  of  correlated  sampling  to  APAIR,  reanalysis 
of  histories  using  weight  factors  to  correct  for  the  difference  in  problem  char¬ 
acteristics  was  used  rattier  than  controlling  the  random  numbers.  Effectively 
the  following  was  performed  for  history  j : 

•  A  base  history  was  run  using  crude  sampling  for  the  long  MAD 
gear  up  to  the  point  where  the  possibility  of  detection  was  to  be 
tested. 

•  Given  the  range,  a  probability  of  detection  for  the  long  range 
gear  (pL)  and  the  short  range  gear  (ps)  were  determined 
from  the  curves  shown  in  Fig.  3. 3. 

•  A  random  number  was  compared  to  pr  to  determine  if  a 
detection  occurred  and  the  history  corifinued,  using  the  results 
of  that  random  decision. 

•  Weighting  factors  were  assigned  to  the  history  according  to 
the  following  rules; 

-  If  the  first  test  for  a  detection  with  the  long  gear  resulted  in 
a  detection,  a  weight  correction  factor  of  -  =  Po/Pt  was 
assigned  to  account  for  the  short  range  gear.  If  me  first 
attempt  resulted  in  no  detection,  he  the  assigned  weighting 
factor  was 


W 


1-PS 


The  history  continues  for  subsequent  tests  of  detection  with  the 
long  MAD  gear  by  assigning  weights  according  to 
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Probability  of  detection 


1.0 


0.5  h 


0 


Long  range  MAD 


Short  range 
MAD 


0.5 


Importance  function 


2.0 


Range  of  Submarine  from  Aircraft 


Fig.  3. 3  Probability  of  detection  versus  range  for  the 
MAD  detectors  used  in  the  demonstration  of 
history  reanalysis  and  importance  sampling. 
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if  a  detection  occurred  on  the  i+1  st  test  and 


1-Po 

w.  -  . « e  •  w.  . 

1+1,  j  1-PL  1,3 

when  a  nondetection  occurred  on  the  i+ 1  st  use  of  the  MAD 
gear.  This  is  continued  until  the  game  stops  due  to  a  kill  or 
some  other  reason  (e.g. ,  sonobuoy  stores  exhausted).  The 
final  weighting  factor  is  defined  as  W^. 

•  If  a  kill  occurred  in  this  history  then 

n.  =  1 
3 

otherwise 

n.  =  0 
1 

(Note  that  the  total  number  of  kills  in  N  histories  is  simply 

vlv 

j=l 

The  above  series  of  steps  was  performed  for  the  N  histories  and  the 
following  formulas  were  used  to  estimate  the  differences  in  the  parameters 
of  interest. 

Consider  first  the  probability  of  kill.  The  estimated  probability  of 
kill  for  the  long  gear  is  given  by 

Pl  '  1/N  X)  "j 

i=i 

The  estimated  probability  of  kill  for  the  short  gear  is  given  by 

N 

is  »  l/N  £  n.W. 

j=l 
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The  difference  in  probability  of  kill  (long  MAD  -  short  MAD)  is 
N 

A  =  1/NVn.  (1-W  ) 

3=1 

Similar  considerations  were  used  for  the  time  to  kill.  If  history  j 
resulted  in  a  kill  and  the  time  of  kill  was  T^,  the  average  time  to  kill  for  the 
long  MAD  is 


n.T. 
J  3 


while  for  the  short  MAD  it  is 


N 


PSN  3=1 


n.T.W. 
3  3  3 


(£sn  is  the  'number*  of  kills  in  the  short  MAD  problem  and  is  the  correct 
normalizing  factor  for  the  time  to  kill  and  other  parameters. )  In  a  similar 
manner  the  remaining  parameters  (number  of  torpedos  or  sonobuoys  used 
per  kill)  were  calculated. 


It  is  clear  the  results  obtained  for  long  and  short  range  detectors  are 
highly  correlated.  The  histories  for  the  two  cases  are  not  only  identical  up 
to  the  point  of  possible  detection,  but  they  continue  to  be  correlated  throughout. 
When  a  kill  is  registered  due  to  the  use  of  the  long  range  gear,  it  is  also 
registered  for  the  short  range  gear  but  carries  a  different  weight.  As  the 
random  choices  made  are  appropriate  to  the  long  gear,  they  will  not  be 
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optimum  for  the  short  range  gear.  This  will  increase  the  variance  of  the 
short  range  gear  estimates,  but  the  high  degree  of  correlation  should  still 
reduce  the  total  variance  of  the  difference. 

A  second  significant  advantage  of  the  above  approach  is  that  the  results 
of  two  problems  have  been  obtained  by  actually  performing  only  one  simulation 
(i.  e. ,  for  the  long  gear)  although  some  additional  computation  is  required  to 
perform  the  weighting.  However,  the  total  computing  time  required  was  only 
53%  of  that  needed  to  perform  two  independent  calculations. 

Thus,  to  generate  the  required  correlated  cases,  approximately  one- 
half  of  the  computational  effort  was  required.  In  addition  to  this  time  saving, 
there  were  substantial  improvements  in  the  variances  of  the  different  parame¬ 
ters  estimated.  These  results  are  summarized  in  Table  3. 11  where  it  is  seen 
that  substantial  improvements  in  the  efficiency  were  realized.  For  example, 
in  estimating  the  stores  used  (torpedos  and  Type  B  buoys)  about  a  factor  of 
10  improvement  in  efficiency  was  obtained.  The  overall  average  efficiency 
was  found  to  be  almost  7.  This  can,  of  course,  be  interpreted  as  a  reduction 
in  computational  time  that  can  be  realized  using  the  correlated  sampling  scheme 
described  here.  ° 

The  'short  range  MAD  detection  probability  curve  used  in  this  technique 
and  that  described  in  the  following  section  is  shown  in  Fig.  3. 3.  This  differs 
from  the  short  range  gear  shown  in  Fig.  3. 2  which  was  used  in  the  other  studies 
of  this  report.  This  change  was  made  because  use  of  the  curve  in  Fig.  3. 2 
would  have  resulted  in  a  considerable  number  of  identically  zero  weights  when¬ 
ever  a  detection  occurred  beyond  the  range  of  the  short  gear.  This  would 
have  destroyed  much  of  the  correlation  in  the  histories  and  counteracted  the 
effect  we  were  trying  to  illustrate.  Therefore,  a  dummy  short  MAD  gear  curve 
having  the  same  limits  of  range  as  the  long  range  gear  but  with  a  lowered  proba¬ 
bility  of  detection  at  the  longer  ranges  was  used.  Paradoxically,  this  had  the 
effect  of  making  the  two  types  of  detector  more  similar,  thus  making  it  harder 
to  calculate  the  difference  between  them. 
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TABLE  3. 11 


Variance  Reduction  Using  History  Reanalysis  For  Estimating  Differences 
Between  Short  and  Long  Range  MAD  Gear 


Estimated  Difference 
in  Expected  Value 
(Long  Range  -  Short 
Range) 

Variance  With 
Crude  Sampling 

Variance  With 
History  Reanalysis 

Efficiency 

Probability  of 
Submarine  Kill 

4. 2  x  10"3 

6. 6  x  10-3 

1.21 

Time  to 

Submarine  Kill 

7.8 

2.05 

7.2 

Number  of  Torpedos 
per  Kill 

1.25  x  10‘2 

2.56  x  10 "3 

9.3 

Sonobuoys  Used  per 
Kill  -  Type  A 

0.097 

0.0289 

6.4 

Sonobuoys  Usea  per 
Kill  -  Type  B 

0.204 

0.040 

9.7 

Average  Efficiency  Im 

provement  z  6. 8 

Ratio  of  time  decrease  with  variance  reduction  arO.53 

History  Reanalysis:  One  run  was  made  using  the  long  range  gear  probability 
distribution  for  detection.  Simultaneously,  calculations  were  made  for  short 
range  gear  using  weighting  factors  to  correct  for  difference  in  probability 
between  the  two  distributions. 


Detection 

<PL> 

No 

Detection 

d-pL) 

Weight 
correction 
for  Short 
Gear 

W-!§ 

PL 

1_ps 

w*irr 

1  ^L 
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The  estimated  differences  in  the  various  parameters  and  their  variances 
were  calculated  by  the  same  batching  techniques  as  described  in  the  previous 
section.  For  this  technique  the  uncorrelated  difference  variances  were  de¬ 
termined  by 

2  a  a  2  a  *2  a  a  2  a 

a  (A)  =  tr  (pL)  +  a  (pg)  »  2  a  (pL) 

where  (pT )  was  determined  from  the  batch  values  of  pT  .  Due  to  the  weight 

corrections,  a  o  (pQ)  calculated  from  the  reanalyzed  histories  would  have 

o  2  A 

been  much  larger  than  or  (pc)  from  an  uncorrelated  case.  Rather  than  make 
a  separate  uncorrelated  run,  it  was  simply  assumed  that  a  (pg)  “cr  (pL)  . 

A  third  variation  involving  correlated  sampling  will  be  described  next. 

3.  8  IMPORTANCE  SAMPLING  WITH  CORRELATION  ESTIMATING 
DIFFERENCES  IN  MAD  GEAR  EFFECTIVENESS 

To  illustrate  the  technique  of  importance  sampling,  which  has  wide 
applicability  at  many  stages  of  a  Monte  Carlo  simulation,  a  demonstration 
involving  an  extension  of  the  previous  correlation  problem  was  devised.  It 
was  very  similar  to  the  calculation  of  the  previous  section  except  that,  in 
place  of  the  long  range  MAD  detection  probability,  an  ’importance  function* 
detection  probability,  mid-way  between  the  long  aid  short  range  curves  as 
shown  in  Fig.  3. 3,  was  used  in  making  detect/no  detect  decisions  in  the  game. 

As  explained  above,  the  use  of  the  long  MAD  curve  to  generate  the  his¬ 
tories  results  in  choices  which  are  not  appropriate  to  the  short  MAD  gear. 

This  increases  the  variance  of  the  short  MAD  part  of  the  calculation.  The  use 
of  an  'importance  function'  curve  generates  choices  which  are  more  optimum 
for  the  short  gear  but  less  optimum  for  the  long  MAD.  It  was  hoped  that 
the  result  would  be  a  reduced  variance  overall. 


The  procedure  used  in  this  simulation  was  a  mere  extension  of  the 
procedure  described  in  the  previous  section.  That  is,  for  history  j  the  follow¬ 
ing  sequence  of  steps  took  place: 

•  A  base  history  was  run  using  crude  sampling.  At  points  in  the 
problem  where  detection  was  to  be  tested,  the  importance  func¬ 
tion  was  used  in  a  random  determination  of  whether  or  not  a 
detection  occurred. 

•  From  the  range,  detection  probabilities  were  derived  from  the 
curves  for  the  'importance'  (pj),  the  long  MAD  (p^),  and  the 
short  MAD  (pg)  as  given  in  Fig.  3. 3. 

•  A  random  number  was  compared  to  Pj  to  determine  if  detection 
occurred. 

•  If  this  first  test  resulted  in  a  detection,  an  assigned  weight 
correction  for  the  short  range  gear  was  set  to 


and  a  weight  for  the  long  range  gear  of 


was  assigned. 

If  no  detection  occurred,  the  respective  weights  assigned  were 


•i-ps 

i-pr 


and 


1-Pi 
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•  The  history  continued  for  subsequent  tests  of  detection  by  assign¬ 
ing  respective  weights  according  to: 


for  detection  on  the  i+lst  test  for  submarine  detection  and 


1_PS 

i-Pi 


for  nondet">ctions  on  the  i±l  st  test  with  similar  equations  for  W^. 
This  is  continued  until  a  kill  occurs  or  the  history  is  otherwise 
terminated.  The  final  weighting  factors  calculated  in  the  history 

are  defined  as  W.  and  W.  . 

]  ] 

•  If  a  kill  occurred  in  this  history,  then 

n.  =  1 
3 

otherwise 

n.  =  0 
3 

(the  total  number  of  kills  in  N  histories  is 


N 


The  above  series  of  steps  was  performed  for  the  N  histories  and  the 
following  formulas  were  used  to  estimate  the  differences  in  the  parameters  of 
interest. 

For  the  probability  of  kill,  the  estimate  for  the  short  range  gear  was 
given  by 


A 


n.WS 
3  3 
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and  for  the  long  range  gear, 


At 


n.WL 
3  3 


The  difference  is  simply 


N 


A  =  PL-Pg  =  1/N^  IjfWj-Wp 

3=1 


Similarly,  if  T.  is  the  time  to  kill  on  history  j  (assuming  a  kill 
occurred),  the  average  time  to  kill  a  submarine  would  be  given  by 

N 


—  n.T.W1 
•  -vr  3  3  3 
pLN  3=1 


L 

3 


and 


N 

p  N 

O  J-X 


n.T.W? 
3  3  3 


The  number  of  torpeaos  ana  sonobuoys  used  were  calculated  in  a  similar 
manner. 


The  correlation  between  the  results  for  the  long  and  short  MAD  detec¬ 
tors  in  this  case  arises  from  the  fact  that  their  estimates  are  derived  from  the 
same  importance  sampled  set  of  histories.  Also,  it  is  important  to  recognize 
that  the  running  time  is  approximately  one -half  of  that  required  to  run  two 
independent  cases.  In  fact,  it  was  found  that  the  ratio  of  the  running  time 
with  and  without  the  correlation  was  a  factor  of  0. 53,  which  is  the  same  as 
the  result  in  the  previous  case. 
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As  expected,  there  was  substantial  improvement  in  the  variance  of  the 
estimates.  The  results  are  summarized  in  Table  3. 12.  It  can  be  seen  that  a 
factor  of  almost  20  was  achieved  in  the  efficiency  of  the  estimator  for  the  dif¬ 
ference  in  the  number  of  torpedos  used.  Also,  a  substantial  improvement  (a 
factor  of  8. 5)  in  the  variance  of  the  difference  in  the  kill  probability  was 
realized.  An  average  efficiency  of  7. 2  was  found  which,  as  before,  can  be  con¬ 
strued  as  a  direct  factor  for  improvement  in  problem  running  time. 
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TABLE  3. 12 


Variance  Reduction  Using  Importance  Sampling  For  Estimating  Differences 

Between  Short  and  Long  Range  MAD  Gear 


Estimated  Differences 
in  Expected  Values 
(Long  Range  -  Short 
Range) 

Variance  With 
Straightforward 
Sampling 

Variance  With 
Correlation 

Efficiency 

Probability  of 
Submarine  Kill 

4.2  x  10  "3 

1. 22  x  10~3 

8.5 

Time  to 

Submarine  Kill 

7.8 

8.3 

1.8 

Number  of  Torpedos 
per  Kill 

1.  2,5  x  10"2 

1.21  x  10  "3 

19.7 

Sonobuoys  Used  per 
Kill  -  Type  A 

0. 097 

0.062 

3.0 

Sonobuoys  Used  per 
Kill  -  Type  B 

0.204 

0.077 

5.0 

Ratio  of  time  decrease  with  variance  reduction  =?0.53 
Importance  Sampling:  One  run  was  made  using  the  importance  function  for 
selection.  Simultaneously  calculations  were  made  for  si,,'-c  and  long  range 
gear  using  weighting  factors  to  correct  for  differences  in  the  probability 
distributions. 


Detection 

<pi> 

No 

Detection 

(i-px) 

Weight 
correction 
(short  gear) 

WS  ps 

*  5 

ws  1_Ps 
w  • 

Weight 
correction 
(long  gear) 

I 

wL  1"Pl 

w  *  I-p~ 
I 
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APPENDIX  A 


MIRAN  -  A  MACHINE  INDEPENDENT  PACKAGE  FOR  GENERATING 
UNIFORM  RANDOM  NUMBERS 


A*  1  GENERAL  DISCUSSION 

The  standard  technique  for  producing  uniform  random  numbers  on 
modern  high-speed  computers  is  an  algorithm  known  as  the  multiplicative 
congruenlial  method.  This  -method  is  expressed  mathematically  as 

Rji+j  =  X.Rn  (modulo  P)  . 


Since  the  R's  are  integers  ranging  from  1  to  P-1,  successive  real  random 

numbers  uniformly  distributed  from  0  to  l  are  generated  by  dividing  Rr  by  P. 

The  properties  of  this  technique  as  a  random  number  generator  (RNG)  are 

highly  dependent  on  the  choice  of  the  generator,  X,  and  the  modulus,  P. 

Unfortunately,  there  are  many  RNGs  in  current  use  which  do  not  approximate 

randomness  closely  enough  to  be  sufficient  for  all  Monte  Carlo  calculations 

and,  what  is  far  worse,  do  manage  to  pass  some  of  the  simple  tests  for 

randomness.  There  are,  however,  several  choices  of  X  and  P  which  have 

(9) 

been  thoroughly  tested,  both  theoretically  and  through  many  years  of  actual 
use  in  Monte  Carlo  calculations,  and  which  appear  to  be  sufficiently  random 
for  general  usage. 


For  reasons  of  convenience  and  efficiency,  P  is  generally  taken  to 
be  2m  where  m  is  the  number  of  bits,  excluding  the  sign  bit,  in  a  single 
word  on  the  particular  computer  being  used.  The  generation  process  starts 
with  a  fixed  generator,  X,  and  a  starting  value,  Rq.  The  full  product 
from  the  multiplication  of  X  and  Rq  would  usually  fill  two  computer  words; 
however,  the  modulo  P  in  the  algorithm  means  that  we  only  need  the  single¬ 
word,  R^,  comprising  the  low  order  half  of  the  X'RQ  product.  The  randua 
number  generation  is  completed  by  converting  R^  to  a  real  variable  and 
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dividing  by  P.  Rj  replaces  RQ  in  storage  in  the  random  number  subroutine 
and  the  process  is  ready  to  begin  anew. 

In  this  sort  of  a  process  there  have  been  two  barriers  to  developing 
a  Fortran  RNG  subroutine  which  would  be  independent  of  the  particular  com¬ 
puter  for  which  it  was  designed.  The  first  is  the  modulus  P,  which  varies 
from  computer  to  computer  as  the  word  length  varies.  [Choosing  a  universal 
value  of  P  to  fit  the  smallest  computer  is  not  a  good  solution  as  the  proper¬ 
ties  of  a  RNG  become  less  random  as  P  is  made  smaller,  to  the  extent  that 
Coveyou  and  MacPherson^  consider  them  questionable  for  P  =  2*** 

OK 

(IBM  360  series)  and  borderline  for  P  =  2  (IBM  7090,  Univac  1108,  etc.).  ] 
The  second  problem  is  that  the  f’gn  bit  of  Rj  may  need  to  be  cleared  follow¬ 
ing  the  multiplication.  Clearing  the  sign  bit  generally  requires  some  trickery 
in  Fortran  which  varies  from  computer  to  computer  as  the  mode  of  represen¬ 
tation  (one's  complement,  two's  complement,  uncomplemented,  etc.)  of 
negative  numbers  varies. 

The  way  around  these  obstacles  is  to  use  an  explicit  multiple  pre¬ 
cision  representation.  The  integers  and  operations  involved  in  the  RNG 
algorithm  are  separated  into  component  pavts  in  such  a  way  that  all  operations 
are  kept  within  a  single  computer  word  and  no  overflows  into  the  sign  bit  are 
made,  thus  avoiding  the  sign-clearing  problem.  Through  multiple  precision 
a  sufficiently  large  modulus  for  good  RNG  properties  may  be  used  even 
though  the  actual  computer  word  size  is  small.  An  initialization  call  must 
be  made  to  convey  to  the  RNG  the  maximum  integer  allowed  on  the  particular 
computer  being  used  so  that  it  can  set  up  an  appropriate  multiple  precision 
representation. 

The  advantage  of  a  RNG  that  is  machine  independent  is  simple:  it 
greatly  facilitates  the  exchange  and  checkout  of  Monte  Carlo  programs  between 
different  computers.  The  price  paid  for  this  advantage  is  also  simple:  it 
is  a  much  slower  method  of  producing  random  numbers.  However,  it  is 
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still  fast  enough  (several  thousand  random  numbers  generated  in  one  second)  that 
the  time  difference  will  not  be  noticed  in  most  Monte  Carlo  applications. 


A.  2  CHOICE  OF  A  SPECIFIC  ALGORITHM  FOR  MIRAN 

The  work  of  Coveyou  and  MacPherson^  has  provided  a  thorough 

theoreMcnl  analysis  of  many  commonly  used  RNGs.  They  show  that  the  cor- 

relaiion  properties  of  a  RNG  are  strongly  dependent  on  the  modulus  P. 

31  35 

For  values  of  P  =  2  or  2  ,  there  must  necessarily  be  a  waviness  or 

graininess  to  the  joint  distribution  of  two,  three,  and  four  consecutive  ran¬ 
dom  numbers  that  could  lead  to  incorrect  results  for  some  Monte  Carlo  cal- 

47 

culations.  For  P  =  2  ,  the  departures  from  true  randomness  are  small 

enough  as  to  be  negligible  for  practical  calculations.  Among  the  specific 
generators,  X,  tested  by  Coveyou  and  MacPherscn,  there  is  one,  X  =  5J 
which  has  good  statistical  properties  and  which  may  be  easily  produced  by 
a  machine  independent  subroutine.  (In  a  subroutine  designed  for  use  on  com¬ 
puters  of  varying  word  length,  specifying  a  fixed  47-bit  integer  through 

15 

data  statements  would  be  difficult.  However,  5  may  easily  be  produced 

by  multiplying  5’s  after  the  exact  multiple  precision  representation  needed 

47  15 

has  been  established.)  In  addition  the  choice  of  P  =  2  and  X=5  has 
an  added  advantage:  this  particular  choice  of  a  RNG  has  seen  long  usage 
(several  thousand  hours  on  a  CDC  1G04  at  Oak  Ridge  National  Laboratory) 
in  Monte  Carlo  computations  without  any  apparent  problems. 
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A.  3  MULTIPLE  PRECISION  REPRESENTATION 

In  the  basic  algorithm  used  by  MIRAN,  X  and  the  R^  values  will 
be  47 -bit  integers.  This  may  exceed  machine  capacity.  To  keep  all  arith¬ 
metic  operations  from  overflowing  a  single  machine  word,  these  integers 
are  stored  in  an  array  wherein  each  word  of  the  array  constitutes  a  'digit1 
in  a  representation  of  the  integer  to  a  particular  base,  This  basis,  called 
BASE,  is  chosen  at  execution  time  so  that  (BASE)  does  not  exceed  the  maxi¬ 
mum  integer  allowed  on  the  particular  computer  being  used.  Thus,  for 
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17 

example,  on  a  machine  with  35 -bit  words  (unsigned),  BASE  would  be  2 
and  each  47-bit  integer  would  be  broken  down  into  3  words  as  follows: 

47 -bit  Integer  Multiple  Precision  Representation 

blb2 . b13^14 - b30b31 - b47  +0 . 0  bl - bl3  word  3 

+  0 . 0  b.,.. .  ..b0_  word  2 

14  30 

+  0 . 0  b^. . . .  b4?  word  1 

Note  that  the  'digits'  are  stored  in  the  array  in  'reverse'  order,  i.e. , 
word  1  is  the  least  significant  17  bits  of  the  number.  Also,  since  17  does 
not  go  evenly  into  47,  the  last  word  contains  only  13  bits. 

Arithmetic  in  a  multiple  precision  representation  is  carried  out  in 
the  same  manner  as  arithmetic  *s  normally  done  by  hand.  The  addition  of 
two  numbers,  for  example,  is  done  digit  by  digit.  When  two  'digits’,  or  words, 
are  added  there  may  be  an  overflow  into  the  18*b  bit  of  the  result.  This  must 
be  detected,  the  overflow  cleared  out,  and  a  carry  of  1  addec  into  the  next 
higher  'digit'.  Multiplication  is  slightly  more  complex.  It  is  again  carried 
out  digit  by  digit  and  the  resulting  products  are  added,  keeping  them  in  appro¬ 
priate  columns,  to  get  the  final  product.  The  multiplication  of  two  ’digits’ 
produces,  of  course,  a  two-digit  product  which  is  initially  contained  in  a 
single  computer  word.  This  must  be  broken  down  into  a  high-order  digit  and 
a  low-order  digit  with  the  high-order  digit  being  added  into  the  next  higher 
column  of  the  result.  As  each  column  is  added,  a  carry  over  into  the  next 
higher  column  may  be  needed.  Thus,  in  our  example  where  three  words  were 
used  for  each  integer,  nine  multiplies  and  several  additions  would  be  needed 
to  form  the  six-word  full  product  as  schematized  below. 


h21  l21 

h31  ^31 

h12  '12 

h22  *22 

h32  l22 

hl3  ^13 

h23  ^23 

h33  ^33 _ _ 

S6  S5  S4  S3  S2  S1 

where  h..  and  l..  are  the  high  and  low  order  parts  of  the  product  of 
U  13 

d.  and  d!  . 

1  3 

A.  4  USE  OF  MIRAN  PACKAGE 
Initialization.! 

Before  generating  any  random  numbers,  it  is  necessary  to  make  an 
initialization  call.  This  is  done  by  the  statement 

CALL  RANSET  (MAXINT,  NSTART) 

where  MAXINT  is  the  maximum  integer  allowed  on  the  computer  (or  compiler) 
being  used.  NSTART  is  the  starting  value,  Rq,  to  be  used  in  the  random 
number  sequence.  If  NSTART  is  less  than  or  equal  to  0 ,  a  default  value 
of  2001  is  supplied  for  NSTART.  If  NSTART  is  even,  the  next  higher  odd 
number  will  be  used. 
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For  example  MAXINT  =  235-l  on  a  1108,  248-l  on  a  CDC-6600,  etc. 
Good  values  for  NSTART  are  any  odd  integer  although  frequent  use  of 
small  odd  integers  is  not  recommended  for  calculations  employing  a  re¬ 
latively  small  number  of  random  numbers. 

The  random  numbers  are  generated  in  subroutine  URAND  which  may 
be  used  as  either  a  function  subroutine  or  as  an  ordinary  subroutine  return¬ 
ing  a  value.  Thus,  either 

CALL  URAND  (R) 
or 

R  =  URAND  (X) 

will  store  a  uniform  random  number  in  R.  (Note  that  in  the  second  form 
the  same  random  number  will  also  be  stored  in  X.  Thus,  X  must  be  a 
Fortran  variable  and  not  a  constant. ) 

Limitations  of  MIRAN: 

MIRAN  will  work  on  all  computers  where  MAXINT  is  greater  than 
94 

1023  and  less  than  2  .  (These  limits  are  practical  and  not  theoretical  and 

could  be  extended  if  it  were  ever  necessary. ) 

A.  5  MIRAN  PROGRAM  DETAILS 

The  Fortran  listings  of  the  two  MIRAN  routines  URAND  and  EANSE-T 
are  presented  in  Figures  A-l  and  A-2.  The  accompanying  logic  flow  is  de¬ 
tailed  in  Figures  A-3  and  A-4.  Additional  explanation  of  the  last  step  in  the 
URAND  logic  is  provided  below. 

The  two  subroutines  URAND  and  RANSET  communicate  through  a 
labelled  common,  MIRNG  which  contains 

RAN(10)  -  An  array  containing  the  'digits'  of  the  current  (or  last) 

multiple  precision  random  integer 
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Figure  A-l.  Fortran  listing  of  URAND 
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Figure  A -2.  Fortran  listing  of  RANSET 
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START 


SND 


Figure  A -3,  Logic  flow  chart  for  URAND 


Figure  A -4.  Logic  flow  chart  for  RANSET 
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GEN(IO)  -  An  array  containing  the  generator  A(=  5  )  in  multiple 
precision  representation 

NWRD  -  The  number  of  words  used  in  the  multiple  precision 
representation  of  an  integer 

BASE  -  The  base  used  in  the  multiple  precision  representation 

MOD  -  The  maximum  value  of  the  highest  order  ’digit*  in  the 

multiple  precision  representation 
FBASE  -  Floating  point  value  of  BASE 

FMOD  -  Floating  point  value  of  MOD 

RAN,  GEN,  NWRD,  and  NBASE  are  Fortran  integers;  FBASE  and  FMOD  are 
Fortran  real  quantities. 

An  alternative  method  (unfortunately,  not  machine  independent)  of  giving 
the  routine  a  starting  value  is  to  save  the  array  RAN  at  the  end  of  a  run  and  to 
restore  RAN  at  the  start  of  the  new  run  (just  after  the  RANSET  call). 

In  the  last  step  of  the  URAND  flow  the  objective  is  conversion 
of  the  multiple  precision  integer  random  number  R  to  a  floating  point 
random  number  X  between  0  and  1.  The  multiplo  precision  integer 
produced  by  the  random  number  algorithm  is  represented  by  the  ’digits’ 
rl’  r2*  * " '  * ,  rn  (remem^ei’  that  r^  is  the  lowest  order  digit.  Thus, 

R  =  rx  +  (BASE)*  rg  +  (BASE)2*r3  + _ +  (BASE)N~*.  r^  . 

Notice  that  we  have,  from  the  manner  in  which  N  and  MOD  were  established. 


P  =  (BASE)**-1.  MOD  . 


The  uniform  random  number  desired  is  given  by  R/P.  Thus  we  have, 


R  _  rl  ,  r2  r3 

*  (BASE)N"1-MOD  (BASE)N“2-MOD  (base)n  -mod 


rN-l  rN 

+ - +  BASE' MOD  +  MOD 


+  .  •  «  • 


BASE  'rl^ - ^ 


MOD  V*N  T  BASE  '*N-1 . BASE  2 

Starting  from  the  right  it  is  easy  to  compute  this  iteratively. 

A.  6  FIRST  100  RANDOM  NUMBERS  PRODUCED  BY  MIRAN 


For  checkout  purposes,  Table  A-l  lists  the  first  100  random  num¬ 
bers  produced  by  MIRAN  when  the  default  value  of  NSTART.,  2001,  is  used 
as  the  starting  random  number. 


First  100  Random  Numbers  Produced  by  Machine -Independent  Random 

Number  Generator 
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